import os
import numpy as np
import pandas as pd
import nibabel as nib

from scipy.spatial.distance import cdist
from sklearn.metrics import pairwise_distances as pdist

import networkx as nx
import plotly.graph_objs as go
import plotly.figure_factory as ff
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import linecache 
from itertools import compress
import xml.etree.ElementTree as ET

import matplotlib as mpl
from matplotlib import cm

# Color maps
viridis = cm.get_cmap('viridis')
tab10   = cm.get_cmap('tab10')

Plot using plotly

## To save as vector image use Orca, e.g.:
# /home/rhaast/.local/bin/orca graph vascular_tree_course.json -o vascular_tree_course -f pdf

Some global plotly settings

# Adjust camera position
camera_vasculature = dict(
    eye=dict(x=0, y=0, z=1.5)
)

camera_hippocampus = dict(
    eye=dict(x=0, y=0, z=2.5)
)

# Hide axis planes
scene = dict(
    xaxis = dict(showbackground=False, visible=False),
    yaxis = dict(showbackground=False, visible=False),
    zaxis = dict(showbackground=False, visible=False)
)

# Set margins
margin = dict(l=0, r=0, b=0, t=0)

# Combine into master dictionary
vasculature_settings = dict(
    width=675,
    height=500,
    scene=scene,
    scene_aspectmode='data',
    scene_camera=camera_vasculature,
    margin=margin,
    showlegend=False,
    hovermode=False
)

hippocampus_settings = dict(
    width=675,
    height=500,
    scene=scene,
    scene_aspectmode='data',
    scene_camera=camera_hippocampus,
    margin=margin,
    showlegend=False,
    hovermode=False
)

IO paths

subject    = os.environ["SUBJECT"]
hemisphere = os.environ["HEMISPHERE"]
# subject    = '01'
# hemisphere = 'R'

mri_maps   = ['B1map','PWI']

vasculature_path = '../results/vasculature/sub-{}/{}'.format(subject, hemisphere)
in_xml           = '{}/vascular_tree.xml'.format(vasculature_path)
in_centerline    = '{}/vascular_tree_centerline.nii.gz'.format(vasculature_path)
in_xfm           = '../resources/mevislabSTL2Orig.txt'

in_lores_surface = '../results/surface_warps/sub-{}/{}/midthickness.downsampled.native.surf.gii'.format(subject, hemisphere)
in_hires_surface = '../results/surface_warps/sub-{}/{}/midthickness.native.surf.gii'.format(subject, hemisphere)

mri_path = '../results/maps/sub-{}'.format(subject, hemisphere)
in_mri   = '{}/sub-{}_{}_{}.nii.gz'

MeVisLab inverse transform

parameters = linecache.getline(in_xfm, 4) 
parameters = parameters.split(' ')[1:10]

xfm = np.array([ float(p) for p in parameters ]).reshape(3,3)

Vascular system graph

tree = ET.parse(in_xml)
root = tree.getroot()

Hippocampal surface

# Lores surface, for visualization
surf             = nib.load(in_lores_surface)
vertices_lores   = surf.darrays[0].data
faces_lores      = surf.darrays[1].data
n_lores_vertices = len(vertices_lores)

# Hires surface, for output maps
surf             = nib.load(in_hires_surface)
vertices_hires   = surf.darrays[1].data
faces_hires      = surf.darrays[0].data
n_hires_vertices = len(vertices_hires)

Centerline NII

nii = nib.load(in_centerline)
centerline = nii.get_fdata()
xfm_nii = nii.affine
voxel_size = abs(nii.affine[0,0])

centerline_idx    = np.argwhere(centerline)
centerline_coords = [ np.dot(np.hstack((centerline_idx[i],1)),xfm_nii.T)[:3] for i in range(0,len(centerline_idx)) ]
centerline_coords = np.array(centerline_coords)

Load MRI maps

mri_data = np.zeros((centerline.shape[0],centerline.shape[1],centerline.shape[2],len(mri_maps)))

for p,parameter in enumerate(mri_maps):
    in_nii   = in_mri.format(mri_path, subject, parameter, hemisphere)
    nii_data = nib.load(in_nii).get_fdata()
    mri_data[1:-1,1:-1,1:-1,p] = nii_data

Extract skeleton

x = []
y = []
z = []
d = []

for skeleton in root.iter('Skeleton'):
    coords = skeleton.find('_PropertyContainer').find('PropertyValues').find('Vector3').text 

    x.append(float(coords.split(' ')[0]))
    y.append(float(coords.split(' ')[1]))
    z.append(float(coords.split(' ')[2]))

    diameter = skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text

    d.append(float(diameter)*2)
            
data_skeleton       = np.vstack((x,y,z,d))
data_skeleton[:3,:] = np.dot(xfm,data_skeleton[:3,:])

Per edge

edge_diameters = []

for skeleton in root.iter('Skeletons'):
    d = []
    for _skeleton in skeleton.iter('Skeleton'):
        diameter = _skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text
        d.append(float(diameter))
    edge_diameters.append(np.median(d))

Extract main nodes

node_id = 0
df_main_nodes = pd.DataFrame(columns=['x','y','z','diameter'])

for node in root.iter('node'):   
    coords = node.find('_BaseGraphItem').find('_PropertyContainer').find('PropertyValues').find('Vector3').text 

    node_x = float(coords.split(' ')[0])
    node_y = float(coords.split(' ')[1])
    node_z = float(coords.split(' ')[2])
    node_xfm = np.dot(xfm,np.array((node_x,node_y,node_z)))    
    
    diameter = node.find('_BaseGraphItem').find('_PropertyContainer').find('PropertyValues').find('double').text 
    
    df_main_nodes.loc[node_id] = [node_xfm[0], node_xfm[1], node_xfm[2], float(diameter)*2]
    
    node_id += 1

Course resolution network edges

Only contains edges between the main (i.e., starting, root and branching) nodes.

# MeVisLab uses 1-based index, so minus one
pre_node  = [int(edge.find('PredId').text)-1 for edge in root.iter('edge')]
post_node = [int(edge.find('SuccId').text)-1 for edge in root.iter('edge')]
                                                                 
df_main_edges = pd.DataFrame(list(zip(pre_node,post_node)), columns=['pre','post'])

Fine resolution network nodes and edges

Uses skeleton elements from XML data

# Initate pandas dataframe for nodes and edges
node_id = 0
df_nodes = pd.DataFrame(columns=['x','y','z','diameter','b1','pwi'])

edge_id = 0
df_edges = pd.DataFrame(columns=['pre','post','distance'])

# Itererate through edges
for edge in root.iter('edge'):

    # Iterate through elements of skeleton
    ids = None
    per_skeleton = None   
    for skeleton in edge.iter('Skeleton'):
        
        # Parse x,y,z coordinates
        skel_xyz = skeleton.find('_PropertyContainer').find('PropertyValues').find('Vector3').text 
        skel_x   = float(skel_xyz.split(' ')[0])
        skel_y   = float(skel_xyz.split(' ')[1])
        skel_z   = float(skel_xyz.split(' ')[2])
        
        # Transform coordinates
        skel_xfm = np.dot(xfm,np.array((skel_x,skel_y,skel_z)))
        
        # Parse diameter
        diameter = skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text 
        
        # Extract MRI-based data
        parameter_data = []
        for p,parameter in enumerate(mri_maps):
            tmp = np.sum(np.abs(centerline_coords - [skel_xfm[0], skel_xfm[1], skel_xfm[2]]),axis=1)
            idx = centerline_idx[np.argmin(tmp)]
            parameter_data.append(mri_data[idx[0],idx[1],idx[2],p])
        
        df_nodes.loc[node_id] = [skel_xfm[0], skel_xfm[1], skel_xfm[2], float(diameter)*2, parameter_data[0], parameter_data[1]]

        # Increase id
        ids = np.hstack((ids, node_id)) if ids is not None else node_id
        node_id += 1

    # Define edges based on collected id's
    for i in ids[:-1]:
        distance = cdist(
            np.array((df_nodes.x[i],df_nodes.y[i],df_nodes.z[i])).reshape(-1,1),
            np.array((df_nodes.x[i+1],df_nodes.y[i+1],df_nodes.z[i+1])).reshape(-1,1)
        )
        df_edges.loc[edge_id] = [i,i+1,distance[0][0]]
        edge_id += 1 

Dataframes need to be cleaned as there will be duplicates near the main nodes

duplicate_entries = df_nodes[df_nodes.duplicated(['x','y','z','diameter'])]

for i in range(0,len(duplicate_entries)):
    duplicates = df_nodes[(df_nodes.x == duplicate_entries.iloc[i].x) &
                          (df_nodes.y == duplicate_entries.iloc[i].y) &
                          (df_nodes.z == duplicate_entries.iloc[i].z) &
                          (df_nodes.diameter == duplicate_entries.iloc[i].diameter)
                         ]
    for d in duplicates.index.values[1:]:
        df_edges = df_edges.replace(d, duplicates.index.values[0])
        df_nodes = df_nodes.drop(d)

Generate networks

Course resolution

G_nodes = df_main_nodes.index.values

G_course = nx.Graph()
G_course.add_nodes_from(G_nodes)
for n in G_nodes:
    G_course.nodes[n]['position'] = (df_main_nodes.x[n], df_main_nodes.y[n], df_main_nodes.z[n])
    G_course.nodes[n]['diameter'] = (df_main_nodes.diameter[n])
    
for e in np.arange(0,df_main_edges.shape[0]):
    G_course.add_edge(
        int(df_main_edges.iloc[e]['pre']),
        int(df_main_edges.iloc[e]['post'])
    )

nodes = np.unique(df_main_nodes.index.values)

G_course_props_dict = {n: (
    df_main_nodes.x[n], df_main_nodes.y[n], df_main_nodes.z[n], df_main_nodes.diameter[n], G_course.degree[n]
) for n in nodes}

Number of edges per main node ('degree')

node_degree_color_dict = {1 : 'green', 2 : 'black', 3 : 'red', 4 : 'purple', 5 : 'blue'}
node_degree_color = [None] * nx.number_of_nodes(G_course)

for node, degree in G_course.degree:
    node_degree_color[node] = node_degree_color_dict[degree]

Visualize

fig = go.Figure()

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    color='lightgray')])

fig.add_trace(go.Scatter3d(
    x=df_main_nodes.x,
    y=df_main_nodes.y,
    z=df_main_nodes.z,
    mode='markers',
    marker=dict(
        size=3,
        color=node_degree_color, #>'black'
    )))

for i,j in enumerate(G_course.edges()):
    x = np.array((G_course_props_dict[j[0]][0], G_course_props_dict[j[1]][0]))
    y = np.array((G_course_props_dict[j[0]][1], G_course_props_dict[j[1]][1]))
    z = np.array((G_course_props_dict[j[0]][2], G_course_props_dict[j[1]][2]))
    fig.add_trace(
        go.Scatter3d(
            x=x,
            y=y,
            z=z,
            mode='lines',
            line=dict(
                width=3,
                color='black'
            )            
        )
    )

fig.update_layout(dict1=vasculature_settings)

fig.write_json("{}/vascular_tree_course.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_course.html".format(vasculature_path))

fig.show()

Fine resolution

Used cleaned dataframes to define network

G_nodes = df_nodes.index.values

G_fine = nx.Graph()
G_fine.add_nodes_from(G_nodes)
for n in G_nodes:
    G_fine.nodes[n]['position'] = (df_nodes.x[n], df_nodes.y[n], df_nodes.z[n])
    G_fine.nodes[n]['diameter'] = (df_nodes.diameter[n])
    
for e in np.arange(0,df_edges.shape[0]):
    G_fine.add_edge(
        int(df_edges.iloc[e]['pre']),
        int(df_edges.iloc[e]['post']),
        distance=float(df_edges.iloc[e]['distance'])
    )

nodes = np.unique(df_nodes.index.values)

G_fine_props_dict = {n: (
    df_nodes.x[n], df_nodes.y[n], df_nodes.z[n], df_nodes.diameter[n], G_fine.degree[n]
) for n in nodes}

Visualize

fig = go.Figure()

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    color='lightgray')])

fig.add_trace(go.Scatter3d(
    x=df_main_nodes.x,
    y=df_main_nodes.y,
    z=df_main_nodes.z,
    mode='markers',
    marker=dict(
        size=3,
        color=node_degree_color, #>'black'
    )))

for i,j in enumerate(G_fine.edges()):
    x = np.array((G_fine_props_dict[j[0]][0], G_fine_props_dict[j[1]][0]))
    y = np.array((G_fine_props_dict[j[0]][1], G_fine_props_dict[j[1]][1]))
    z = np.array((G_fine_props_dict[j[0]][2], G_fine_props_dict[j[1]][2]))
    fig.add_trace(
        go.Scatter3d(
            x=x,
            y=y,
            z=z,
            mode='lines',
            line=dict(
                width=10,
                color='black'
            )            
        )
    )
    
fig.update_layout(dict1=vasculature_settings)

fig.write_json("{}/vascular_tree_fine.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_fine.html".format(vasculature_path))

fig.show()

Calculate volume and surface area for vessel segments

Edges are treated as truncated cones with varying radii and height (i.e., radii and distance between adjacent nodes)

def calculate_edge_geometry(r1, r2, h):
    import math
    v  = (1/3)*math.pi*(pow(r1,2)+r1*r2+pow(r2,2))*h # volume
    la = math.pi*(r1+r2)*np.sqrt(pow((r1-r2),2)+pow(h,2)) # lateral area
    sa = la + math.pi*(pow(r1,2)+(pow(r2,2))) # surface area
    
    return(v, la, sa)

edge_geometry = [calculate_edge_geometry(
    df_nodes.diameter[df_edges.pre[i]]/2, df_nodes.diameter[df_edges.pre[i]]/2, df_edges.distance[i]
) for i in df_edges.index.values]
edge_geometry = np.array(edge_geometry)

df_edges['volume']       = edge_geometry[:,0]
df_edges['lateral_area'] = edge_geometry[:,1]
df_edges['surface_area'] = edge_geometry[:,2]

# Project edge geometry data on neighboring nodes
node_geometry = np.zeros((nx.number_of_nodes(G_fine),6))

for n, node in enumerate(df_nodes.index.values):
    df_filtered = df_edges[(df_edges.pre==node) | (df_edges.post==node)]
    
    # Vessel geometry
    node_geometry[n,0] = np.mean(df_filtered[df_filtered.volume!=0].volume)
    node_geometry[n,1] = np.mean(df_filtered[df_filtered.volume!=0].lateral_area)
    node_geometry[n,2] = np.mean(df_filtered.surface_area)
    node_geometry[n,3] = np.mean(
        df_nodes.loc[pd.unique(
            df_filtered[['pre','post']].values.ravel()
        )].diameter
    )
    
    # Smooth MRI-sampled data
    node_geometry[n,4] = np.mean(
        df_nodes.loc[pd.unique(
            df_filtered[['pre','post']].values.ravel()
        )].b1
    )
    node_geometry[n,5] = np.mean(
        df_nodes.loc[pd.unique(
            df_filtered[['pre','post']].values.ravel()
        )].pwi
    )      

# Transform to hex codes for plotting
volume_norm  = mpl.colors.Normalize(vmin=min(node_geometry[:,0]), vmax=max(node_geometry[:,0]))
volume_color = [ mpl.colors.to_hex(viridis(volume_norm(node_geometry[i,0]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]

lateral_area_norm  = mpl.colors.Normalize(vmin=min(node_geometry[:,1]), vmax=max(node_geometry[:,1]))
lateral_area_color = [ mpl.colors.to_hex(viridis(lateral_area_norm(node_geometry[i,1]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]

surface_area_norm  = mpl.colors.Normalize(vmin=min(node_geometry[:,2]), vmax=max(node_geometry[:,2]))
surface_area_color = [ mpl.colors.to_hex(viridis(surface_area_norm(node_geometry[i,2]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]

Normalize diameter data for plotting

from sklearn.preprocessing import minmax_scale
diameter_scaled = minmax_scale(node_geometry[:,3], feature_range=(3,10))

Connected components

Identify connected components in graph and sort based on size

# Identify connected components and sort from largest to smallest (i.e., connected nodes)
G_cc = sorted(nx.connected_components(G_fine), key=len, reverse=True)

# Treshold for masking small components
treshold = 200
components_mask = np.ones((nx.number_of_nodes(G_fine)))

# Create compononts color map
components_color = [None] * nx.number_of_nodes(G_fine)
for c, component in enumerate(G_cc):
    G_component     = G_fine.subgraph(G_cc[c])
    intersect, ind_a, ind_b = np.intersect1d(
        list(G_fine.nodes), list(G_component.nodes), return_indices=True)

    for i in ind_a:
        components_color[i] = mpl.colors.to_hex(tab10(c))
        
    # Mask if components smaller than treshold
    if len(ind_a) <= 100:
        components_mask[ind_a] = 0

components_mask = np.array(components_mask,dtype='bool')

Visualize

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    color='lightgray',   
    showscale=False)])

fig.add_trace(go.Scatter3d(
    x=df_nodes.x[components_mask], 
    y=df_nodes.y[components_mask], 
    z=df_nodes.z[components_mask],
    mode='markers',
    marker=dict(
        size=diameter_scaled[components_mask],
        line_width=0,
        opacity=1,
        color=node_geometry[components_mask,3], #list(compress(node_geometry[:,3], components_mask))
        colorscale='plasma',
        showscale=True,
        colorbar=dict(
            title="Diameter (mm)",
            titleside="right",
            thickness=10
            )        
    )))

fig.update_layout(dict1=vasculature_settings)

fig.write_json("{}/vascular_tree_diameter.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_diameter.html".format(vasculature_path))

fig.show()

Shortest path

Generate shortest path from root nodes and export to matrix for color mapping

# Extract root nodes from XML
root_nodes = [int(r.text) for r in root.iter('RootID')]

# Identify root nodes based on x, y, z coordinates 
root_node_ids = []
for root_node in root_nodes:
    x = df_main_nodes.iloc[root_node-1].x
    y = df_main_nodes.iloc[root_node-1].y
    z = df_main_nodes.iloc[root_node-1].z
    
#     x = df_main_nodes[1,df_main_nodes[0,:]==root_node][0]
#     y = df_main_nodes[2,df_main_nodes[0,:]==root_node][0]
#     z = df_main_nodes[3,df_main_nodes[0,:]==root_node][0]
    
    df_root = df_nodes.query('x=={} & y=={} & z=={}'.format(x,y,z))
    
    root_node_ids.append(df_root.index.values[0])

# Initiate array for filling with calculated distances
distance_to_root = np.zeros((len(G_fine.nodes),len(root_node_ids)))

for r, root_node in enumerate(root_node_ids):
    # Define source node
    source_node = root_node

    # Calculate shortest path from source node to each other connected node
    shortest_path_from_root = nx.shortest_path(G_fine, source=source_node, weight='distance')

    # Iterate through each connected node and its path to the source node
    for node, path in shortest_path_from_root.items():

        # Find index of node
        node_index = np.where(G_nodes == node)[0][0]

        # Skip source node
        if node != source_node:

            mask = list(zip(path[:-1], path[1:]))

            df_x = pd.MultiIndex.from_frame(df_edges[['pre','post']])
            df_y = pd.MultiIndex.from_frame(df_edges[['post','pre']])

            dist_x = sum(df_edges[df_x.isin(mask)].distance.values)
            dist_y = sum(df_edges[df_y.isin(mask)].distance.values)

            distance_to_root[node_index,r] = dist_x + dist_y        

Heuristically define starting node of largest component

Most anterior, inferior vessel segment characterized by largest diameter

start_nodes = np.zeros((len(G_cc)))
for c, component in enumerate(G_cc):
    order_x = df_nodes['x'][list(G_cc[c])].argsort()
    rank_x  = order_x.argsort()

    order_y = df_nodes['y'][list(G_cc[c])].argsort()[::-1]
    rank_y  = order_y.argsort()

    order_d = df_nodes['diameter'][list(G_cc[c])].argsort()[::-1]
    rank_d  = order_d.argsort()
    
    # x, y, and diameter are weighted differentely. 
    rank_sum       = np.average(np.vstack((rank_x,rank_y,rank_d)).T,weights=[0.5,1,1], axis=1)
    start_nodes[c] = df_nodes.loc[list(G_cc[c])].iloc[np.argmin(rank_sum)].name

Calculate distance from starting node

distance_to_start = np.zeros((len(G_fine.nodes),len(G_cc)))

for s, start_node in enumerate(start_nodes):
    # Define source node
    source_node = start_node
    
    # Calculate shortest path from source node to each other connected node
    shortest_path_from_root = nx.shortest_path(G_fine, source=source_node, weight='distance')

    # Iterate through each connected node and its path to the source node
    for node, path in shortest_path_from_root.items():

        # Find index of node
        node_index = np.where(G_nodes == node)[0][0]

        # Skip source node
        if node != source_node:

            mask = list(zip(path[:-1], path[1:]))

            df_x = pd.MultiIndex.from_frame(df_edges[['pre','post']])
            df_y = pd.MultiIndex.from_frame(df_edges[['post','pre']])

            dist_x = sum(df_edges[df_x.isin(mask)].distance.values)
            dist_y = sum(df_edges[df_y.isin(mask)].distance.values)

            distance_to_start[node_index,s] = dist_x + dist_y  

Visualize

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    color='lightgray', 
    showscale=False)])

fig.add_trace(go.Scatter3d(
    x=df_nodes.x[components_mask], 
    y=df_nodes.y[components_mask], 
    z=df_nodes.z[components_mask],
    mode='markers',
    marker=dict(
        size=diameter_scaled[components_mask],
        line_width=0,
        opacity=1,
        color=np.max(distance_to_start[components_mask],axis=1),
        colorscale='Viridis',
        cmax=60,
        cmin=0,
        showscale=True,
        colorbar=dict(
            title="Geodesic distance from seed (mm, per component)",
            titleside="right",
            thickness=10
            )
    ),

))

fig.update_layout(
    dict1 = vasculature_settings,
    scene=dict(
        xaxis = dict(showbackground=False, visible=False),
        yaxis = dict(showbackground=False, visible=False),
        zaxis = dict(showbackground=False, visible=False),
        annotations=[
            dict(
                x=df_nodes.loc[start_nodes[0]]['x'],
                y=df_nodes.loc[start_nodes[0]]['y'],
                z=df_nodes.loc[start_nodes[0]]['z'],
                text="Seed",
                bgcolor="black",
                font=dict(
                    color="white",
                    size=16
                ),
                showarrow=True,
                arrowcolor="black"
            )
        ])
)

fig.write_json("{}/vascular_tree_distance_to_seed.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_distance_to_seed.html".format(vasculature_path))

fig.show()

Calculate distance between hippocampus and vasculature

hippocampus_vessel_metrics_hires = np.zeros((n_hires_vertices,3))
hippocampus_vessel_metrics = np.zeros((n_lores_vertices,3))

# Hires for exporting
for i in range(0,n_hires_vertices):
    dist = cdist(vertices_hires[i].reshape(1,3), df_nodes.values[:,:3])
    hippocampus_vessel_metrics_hires[i,0] = df_nodes.values[:,3][np.argmin(dist)]
    hippocampus_vessel_metrics_hires[i,1] = np.max(distance_to_start,axis=1)[np.argmin(dist)]
    hippocampus_vessel_metrics_hires[i,2] = np.amin(dist)-(0.5*hippocampus_vessel_metrics_hires[i,0])
    
# Lores for exporting
for i in range(0,n_lores_vertices):
    dist  = cdist(vertices_lores[i].reshape(1,3), df_nodes.values[:,:3])
    hippocampus_vessel_metrics[i,0] = df_nodes.values[:,3][np.argmin(dist)]
    hippocampus_vessel_metrics[i,1] = np.max(distance_to_start,axis=1)[np.argmin(dist)]
    hippocampus_vessel_metrics[i,2] = np.amin(dist)-(0.5*hippocampus_vessel_metrics[i,0])

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    intensity=hippocampus_vessel_metrics[:,2],
    colorscale='magma',    
    showscale=True,
    colorbar=dict(
        title="Euclidean distance to closest vessel (mm)",
        titleside="right",
        thickness=10
        )
)])

fig.update_layout(dict1=hippocampus_settings)

fig.write_json("{}/hippocampus_distance_to_vessel.json".format(vasculature_path))
fig.write_html("{}/hippocampus_distance_to_vessel.html".format(vasculature_path))

fig.show()

Map MRI data onto hippocampus

hippocampus_mri_metrics_hires = np.zeros((n_hires_vertices,len(mri_maps)))
hippocampus_mri_metrics = np.zeros((n_lores_vertices,len(mri_maps)))

# Hires for exporting
for i in range(0,n_hires_vertices):
    idx = np.round(np.dot(np.hstack((vertices_hires[i,:],1)),np.linalg.inv(xfm_nii).T)) 
    for p,parameter in enumerate(mri_maps):
        hippocampus_mri_metrics_hires[i,p] = mri_data[int(idx[0]),int(idx[1]),int(idx[2]),p]

# Lores for exporting
for i in range(0,n_lores_vertices):
    idx = np.round(np.dot(np.hstack((vertices_lores[i,:],1)),np.linalg.inv(xfm_nii).T)) 
    for p,parameter in enumerate(mri_maps):
        hippocampus_mri_metrics[i,p] = mri_data[int(idx[0]),int(idx[1]),int(idx[2]),p]

Visualize

fig = go.Figure(data=[go.Mesh3d(
    x=vertices_lores[:,0], 
    y=vertices_lores[:,1], 
    z=vertices_lores[:,2], 
    i=faces_lores[:,0], 
    j=faces_lores[:,1], 
    k=faces_lores[:,2], 
    intensity=hippocampus_mri_metrics[:,0],
    colorscale='viridis',
    cmax=6,
    cmin=15,    
    showscale=False)])

fig.add_trace(go.Scatter3d(
    x=df_nodes.x[components_mask], 
    y=df_nodes.y[components_mask], 
    z=df_nodes.z[components_mask],
    mode='markers',
    marker=dict(
        size=diameter_scaled[components_mask],
        line_width=0,
        opacity=1,
        color=node_geometry[components_mask,4],
        colorscale='viridis',
        cmin=6,
        cmax=15,
        showscale=True,
        colorbar=dict(
            title="B1+ (mT)",
            titleside="right",
            thickness=10
        )
    )))

fig.update_layout(dict1=vasculature_settings)

fig.write_json("{}/vascular_tree_b1.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_b1.html".format(vasculature_path))

fig.show()

Correlations among vascular morphology

import seaborn as sns

axis_labels = ['Volume (mm$^3$)','Lateral area (mm$^2$)','Surface area (mm$^2$)']

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10,4))

for g in np.arange(0,3):
    sns.scatterplot(x=node_geometry[:,3], y=node_geometry[:,g], alpha=0.25, linewidths=0, ax=axes[g])
    sns.despine()
    
    axes[g].set_xlabel('Diameter (mm)')
    axes[g].set_ylabel(axis_labels[g])

plt.tight_layout()

Save data to file

# Save network to Panda's pickle
df_nodes.to_pickle('{}/sub-{}_{}_network_nodes.pkl'.format(vasculature_path, subject, hemisphere))
df_edges.to_pickle('{}/sub-{}_{}_network_edges.pkl'.format(vasculature_path, subject, hemisphere))

# Save others to Numpy npy
np.save('{}/sub-{}_{}_network_components.npy'.format(vasculature_path, subject, hemisphere),G_cc)
np.save('{}/sub-{}_{}_hippocampus_vessel.npy'.format(vasculature_path, subject, hemisphere),hippocampus_vessel_metrics_hires)
np.save('{}/sub-{}_{}_hippocampus_mri.npy'.format(vasculature_path, subject, hemisphere),hippocampus_mri_metrics_hires)

And done! 💪